#include "../include/Multigrid_Solver.h"
#include "../include/Multigrid_Solver_1D.h"
#include "../include/Multigrid_Solver_2D.h"
#include "../include/Multigrid_Factory.h"
#include "../include/Functions_test.cpp"
#include <string>

std::string Rtype2str(int i)
{
    if (i == 1)
        return "FullWeight";
    else if (i == 2)
        return "Injection";
    else
        return "unkown";
};

std::string Itype2str(int i)
{
    if (i == 1)
        return "Linear";
    else if (i == 2)
        return "Quadratic";
    else
        return "unkown";
};

std::string Stype2str(int i)
{
    if (i == 1)
        return "WTJacobi";
    else if (i == 2)
        return "GaussSeidal";
    else
        return "unkown";
};

int main()
{
    int maxium_iteration = 100;
    std::cout << "Nonhomogeneous Possion Equation:" << std::endl;
    std::cout << "u''(x) = -(sin(x)-cos(x)*cos(x))*exp(sin(x)) ; u(0) = 1, u(1) = exp(sin(1))." << std::endl;
    int numofgrid;
    double _a[1] = {0};

    VectorXd solution_temp, rhs, r_temp, r_origin, solution_exact, sol_err_box, res_err_box;
    double r_origin_norm;
    Factory1DMGSolver Solver1D;
    FactoryMGSolverBase *_p_fac_base = &Solver1D;
    Multigrid_Solver *_arena1D = _p_fac_base->CreateMGSolver();
    int flag1 = 0, flag2 = -1;
    MatrixXd _Convergence_Rate_Save;
    _Convergence_Rate_Save.setZero(8, 4);
    numofgrid = 128;
    int cnt;
    double h;
    char *Rstr, Istr;
    flag1 = 0, flag2 = -1;
    for (numofgrid = 128; numofgrid <= 1024; numofgrid *= 2)
    {
        flag1 = 0;
        flag2 += 1;
        for (int _Stype = 1; _Stype <= 2; _Stype++)
        {
            for (int _Rtype = 1; _Rtype <= 2; _Rtype++)
            {
                for (int _Itype = 1; _Itype <= 2; _Itype++)
                {
                    _arena1D->Initialize(numofgrid, Dim1_type1_df2rd, Dim1_type1_f);
                    _arena1D->setRestircType(_Rtype);
                    _arena1D->setIntpType(_Itype);
                    _arena1D->setSmoothType(_Stype);
                    solution_temp = _arena1D->get_u_origin();
                    rhs = _arena1D->get_f_origin();
                    h = _arena1D->get_h();
                    r_origin = _arena1D->ResidualSolver(solution_temp, rhs, h);
                    r_origin_norm = r_origin.lpNorm<Eigen::Infinity>();
                    solution_exact.setZero(numofgrid + 1);
                    for (int i = 0; i <= numofgrid; i++)
                    {
                        _a[0] = 1.0 * i / numofgrid;
                        solution_exact(i) = Dim1_type1_f(_a);
                    };
                    sol_err_box.setZero(maxium_iteration + 2);
                    res_err_box.setZero(maxium_iteration + 2);
                    sol_err_box(0) = solution_exact.lpNorm<Eigen::Infinity>();
                    res_err_box(0) = r_origin_norm;
                    for (cnt = 1; cnt <= maxium_iteration; cnt++)
                    {
                        solution_temp = _arena1D->Vcycle(solution_temp, rhs, h);
                        r_temp = _arena1D->ResidualSolver(solution_temp, rhs, h);
                        sol_err_box(cnt) = (solution_temp - solution_exact).lpNorm<Eigen::Infinity>();
                        res_err_box(cnt) = (r_temp).lpNorm<Eigen::Infinity>();
                        if (res_err_box(cnt) <= r_origin_norm * 1e-8)
                            break;
                    };
                    _Convergence_Rate_Save(flag1, flag2) = sol_err_box(cnt);
                    flag1++;
                    std::cout << "### 1D_VC:" << numofgrid << "," << _Stype << "," << _Rtype << "," << _Itype << "### (numofgrid,Stype,Rtype,Itype)" << std::endl;
                    std::cout << "This V-cycle with "
                              << "#n = " << numofgrid << " needs #" << cnt << " iterations for 1e-8 re_acu , when :" << std::endl;
                    printf("#Smoother               #Restricter             #Interpolater \n");
                    printf("%10s    %20s    %20s\n", Stype2str(_Stype).c_str(), Rtype2str(_Rtype).c_str(), Itype2str(_Itype).c_str());
                    std::cout << "#Detail:" << std::endl;
                    printf("%17s %15s %15s %15s %15s\n", ((std::string) "Relative_Accuracy").c_str(), ((std::string) "Residual").c_str(), ((std::string) " Ratio_Residual").c_str(), ((std::string) "SolutionError").c_str(), ((std::string) "Ratio_error").c_str());
                    for (int i = 1; i <= cnt; i++)
                    {
                        printf("%2d %14.2e %14.2e %14.2e %14.2e %14.2e\n", i,  res_err_box(i) / res_err_box(0), res_err_box(i),res_err_box(i) / res_err_box(i - 1), sol_err_box(i), sol_err_box(i) / sol_err_box(i - 1));
                    };
                    std::cout << " " << std::endl;
                };
            };
        };
    };

    MatrixXd _Generater_convergence;
    _Generater_convergence.setZero(8, 3);
    for (int j = 0; j <= 2; j++)
    {
        for (int i = 0; i <= 7; i++)
        {
            _Generater_convergence(i, j) = -log2(_Convergence_Rate_Save(i, j + 1) / _Convergence_Rate_Save(i, j));
        };
    };

    std::cout << "##Convergence Level Table" << std::endl;
    MatrixXd A(1, 4);
    A << 0, 128, 256, 512;
    MatrixXd B(9, 1);
    B << 0, 111, 112, 121, 122, 211, 212, 221, 222;
    MatrixXd C(9, 4);
    C.row(0) = A;
    C.col(0) = B;
    for (int j = 1; j <= 3; j++)
        for (int i = 1; i <= 8; i++)
            C(i, j) = _Generater_convergence(i - 1, j - 1);
    std::cout << C << std::endl;
    std::cout << "PS:The first col is the combination code of IntpType,RestrictType and SmootheType." << std::endl;
    std::cout << "Interpolation Function: Linear 　= 　1；Quadratic 　= 　2；" << std::endl;
    std::cout << "Restriction Function:   FullWeighting 　= 　1； Injection 　= 　2；" << std::endl;
    std::cout << "Smoothing Function: WeightedJacobi ,omega = 2/3　= 　1;   GaussSeidel 　= 　2；" << std::endl;
    std::cout<<"EOF"<<std::endl;

    return 0;
};